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Abstract. — We consider the Klein-Gordon equation in the non-relativistic 
limit regime, i.e. the speed of light c tending to infinity. We construct an 
asymptotic expansion for the solution with respect to the small parameter 
depending on the inverse of the square of the speed of light. As the first terms 
of this asymptotic can easily be simulated our approach allows us to construct 
numerical algorithms that are robust with respect to the large parameter c 
producing high oscillations in the exact solution. 



The Klein-Gordon equation is a fundamental physical equation describing 
the motion of a spin-less particle, see for instance [4], [22] and [5] for a physical 
introduction on this topic and [13], [6] and the references therein for the 
mathematical analysis, and in particular [3, 1] for some long time results 
when the equation is set on a compact domain, the case mainly considered in 
the present work. 

In the relativistic regime, i.e. the speed of light c = 1, the Klein-Gordon 
equation (as a non-linear wave equation) is numerically well studied, see for 
instance [23], [16], [20], and, for long time behavior results in the weakly 
nonlinear regime, [7], [10] and [11, 12]. Here we are interested in the Klein- 
Gordon equation in the non-relativistic limit regime, i.e. the speed of light c 
tending to infinity. Thus, our model problem reads 



where we will consider essentially the case f{z) = A|zp^z, A G M, assume 
that the initial condition (p and 7 do not depend on c, and consider periodic 
boundary conditions in the space variable. In this regime the Klein-Gordon 
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equation has recently gained a lot of attention, see for instance [18], [19] and 
[24] (when Equation (1.1) is set on the whole space M'^). Numerically, the 
highly oscillatory nature of the exact solution is very challenging as standard 
numerical schemes require severe time step restrictions depending on the large 
parameter c^, see for instance [14], [8] and [21] for a numerical overview on 
highly oscillatory problems. In the recent work [2], Gautschi-type exponential 
integrators (proposed in [15] for equations with high frequencies generated by 
the linear part) were applied to the non-relativistic Klein-Gordon equation, 
which allow time steps of order 0{c~'^). 

Following the one-term asymptotics results given in [19] we study in this 
paper the existence of a complete asymptotic expansion of the exact solution of 
(1.1) in terms of over a c -independent interval (0, T), and for smooth initial 
data. This will allow us to construct numerical schemes which do not need 
to obey any c-dependent smallness restriction. In a nutshell, the asymptotic 
expansion up to the first-order correction term of (1.1) reads 



(1.2) z{t,x) = luo{t,x)e''''^ + ^vo{t,x)e~''''^ + 0{c-'^), t G [0,T] 



The essential point is that the highly-oscillating part in the asymptotic expan- 
sion (1.2) is only contained in the phases e^^ * and e~^'^ *, but does not appear 
in the above Schrodinger equations for uq and vq. In particular the numerical 
time integration of (1.3) can be carried out without any c-dependent time step 
restriction. Thus, with this asymptotic expansion we are able to obtain a well 
suited numerical approximation to the exact solution z{t) of the Klein-Gordon 
equation in the non-relativistic limit c ^ 1, without any time step restriction 
as we only need to solve the Schrodinger equations (1.3) and multiply the 
numerical approximations of uq and vq with the highly-oscillatory phases e*'^ * 
and e~*'^ *, respectively. 

Note that in comparison with [19] and motivated by numerical implemen- 
tation, we choose to work here with periodic boundary conditions, that is the 
space variable x G T'^ the d-dimensional torus. Moreover, we show that for 
sufficiently smooth initial values the asymptotic expansion of the exact so- 
lution to (1.1) exists up to an arbitrary order, and can be reconstructed by 
superposing highly oscillatory terms to cascaded solutions of c-independent 
Schrodinger-like systems. After a full discretization of these cascaded terms, 
we are able to build numerical schemes which approach the exact solution up 



where uq and vq solve the coupled nonlinear Schrodinger equations 




t;g(0) =ip + i'y. 



Uo{0) =ip-ij, 
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to errors terms of order 0{t^ + c~^^~^ + h^) by simulating the terms in the 
asymptotic expansion up to order with a numerical scheme of order k in 
time (with step size r) and Fourier pseudo-spectral discretization on a grid 
with mesh size h. Note that here, h and r can be chosen independently of c, 
and that s = s{N) depends on the regularity of the solution and of the degree 
of the asymptotic expansion. 

The paper is organized as follows. We commence in Section 2 with the 
essential a priori bounds on the exact solution of the Klein-Gordon equation 
(1.1) showing the existence on a time interval independent of the large pa- 
rameter c. In Section 3 we will derive an asymptotic expansion of the exact 
solution in powers of c~^. In Section 4 we will precisely state the first terms in 
this expansion, which will allow us to construct numerical schemes which are 
robust with respect to c, see Section 6. Furthermore, in Section 5 we analysis 
the asymptotic behavior of the conserved physical quantities. 



2. A priori bounds on the exact solution 

We consider our model problem (1.1) set on the d-dimensional torus T*^. 
We make the following assumptions: 

Assumption 2.1. — The initial values ip and 7 are smooth functions inde- 
pendent of c. Furthermore we assume that the nonlinearity / is gauge invariant 
and polynomial of degree r > 1 in 2; and z. This means that for all number 
a G M and all z £ C, we have /(e*°z) = e*"/(z). Eventually, we assume that 
/ is real in the sense that for all z £ C, we have f{z) = f{z). In particular, 
this shows that if z{0) and dtz{0) are in M, then the solution z{t) is real for 
all times. Typical examples of such nonlinearity are given by the polynomials 
f{z) = A|z|2Pz, p>l. 

For ease of presentation and numerical motivations, we will mainly work on 
the d-dimensional torus T"^, d € N\{0}. However, as the only arguments used 
in the sequel are the local Lipschitz nature of the nonlinearity in Sobolev space 
and some explicit calculations for diagonal operators in Fourier, the reader 
should be convinced that all the result of this paper can be easily extended to 
functions defined on W^. 

For a function u{x), x G T"* in := if*(M,C) we set 

\\u\\i^ = Y.{i+\k\r\uk\', 

where 

n. = ^/^^e-'^-Mx)d., 
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k G Z'^, denote the Fourier coefficients associated with u. Here, we have set 

k ■ X = kixi + • • • + kdXd and |/c|^ = kf + ■ ■ ■ + k'^, 

for k = {ki, . . . , kd) G Z'^ and x = {xi, . . . , Xd) G T'^. For a given c > 0, we 
define the operator 

(V)c = \/-A + c2. 

In Fourier, the operator (V)c is represented by the diagonal operator with 
coefficients 

{{V)c)ki = SkeV\k\^ + c^ K ^ G Z'', 
where 5kt denotes the Kronecker operator. In particular note that the oper- 
ator c(V)~^ is uniformly bounded with respect to c: we have for all k G Z'^, 
\{c(y)^^)kk\ < 1, which implies that for all functions u, ||c(V)^"'^n||^^ < 
With these notations. Equation (1.1) can be written as 

(2.1) dttz + c^{V)lz = c^f{z). 

In order to rewrite the above equation as a first-order system in time, we set 

(2.2) u = z-ic-^{V)^^dtz, v = z - ic-\V)-^dtz. 
Note that we have 

(2.3) z = l(^u + v), 

and if z is real, then u = v. A short calculation shows that in terms of the 
variables u and v equation (2.1) reads 

idtu = -c{V),u + c{V)-'f{l{u + v)), 

(2.4) 

idtv = -c{V)cV + c{V)-'f{^{u + v)). 
The initial boundary conditions associated with this problem are (see (1-1)) 



(2.5) u(0,x) = V5(x) - ic(V);:S(a;), and v{0, x) = (p{x) - ic{V)~'^'y{x). 
Eventually, we set 

(2.6) w=[ , F{w)=\ , and ^= _ . ,^.,_] , 

so that the problem (2.4)-(2.5) can be written as 

(2.7) idtw = -c{V)cW + c{V)^^F{w), w{0) = tp. 

We see that this equation is posed on x H^. For an element w = [u, vY' G 
X H^, we set 
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In the previous and in the following, if an operator - like (V)c - acts on H'^, 
we naturally extend diagonally its action on H'^ x by setting for instance 
for w = {u, v)'^, 

(v,...,v,(:;).p;). 

The solution of the linear equation = in the previous equation) can be 
written as follows: For w = {u,v) G H'^ x , 

w{t) = e^*"<^>-^, 

where the right hand-side is easily defined in Fourier by the formula 

where ipj^ denote the (vector valued) Fourier coefficients of if). Note that the 
linear flow is an isometry in the sense that for ■0 G x H'^, 

Note that in the sequel, and because the equation is set on the torus, we will 
not use any argument involving Strichartz estimates as in [19] but we will 
require more regularity on the initial conditions 7 and (p. 

As we assume that the function / is polynomial of degree r > 2 in z and z 
the application F satisfies estimates of the form 

(2.8) \\F{wi) - F{W2)\\^, < C{1 + \\wi\\^^ + \\w2\\^J-^\\wi - w^W^, 

for s > d/2, with some constants C, and for wi and W2 in x . These 
estimates are consequences of classical bilinear estimates in H'^: 

(2.9) \\uv\\^, <Cs\\u\\^J\v\\^^ for s > d/2, 

for some constants Cs- Moreover the assumption that / is real implies that 
for all we H" X H", F{w) = F{w). 

With these notations, the mild formulation of (2.7) is given by 

w{t) = e**^<^>=^/j + [\{V)-^e'^*~'^^^^'^^F{w{s))ds. 
Jo 

Using the fact that the linear flow e'^^^^^)^ jg isometry in and that c(V)^"'^ 
is a bounded operator in , in particularly uniformly bounded with respect 
to c, we easily derive the following result by a standard fix point argument 
combined with (2.8): 

Proposition 2.2. — Let s > d/2 and M > a given constant. Then there 
exist constants Tg > and Bg > such that for all c > and for all ip G 
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H'^ X such that < M, there exists a mild solution t i— >■ w{t) G 

C^{{0,Ts), X H^) to the equation (2.7), and such that 

vtG(o,r,), \\wmHs<Bs. 

In terms of u and v, the previous proposition can be reformulated as 

Corollary 2.3. — Let s > d/2, and assume that ip G and 7 G H" are 
given. Then there exist constants > and Bs > such that for all c there 
exists a mild solution {u{t),v{t)) € C^{{0,Ts), x H^) to the equation (2.4) 
satisfying the initial boundary conditions (2.5) and such that 

VtG(0,r,), \\u{t)\\^^ + \\v{t)\\^,<Bs. 

Proof. — Using the expression (2.6) for in terms of (p and 7, and the fact 
that c(V)~"^ is uniformly bounded in with respect to c, we obtain that lliAH^^ < 
M, a constant independent of c. This shows that Proposition 2.2 applies. □ 

This corollary immediately gives the existence of the exact solution z(t) of 
(1.1) on a time interval (0,Ts) and with iJ^'-bounds both independent of c. 



3. Asymptotic expansion 

On the c-independent time interval (0, Tg) given by the previous results, we 
search w (at least formally) under the form 

f[J{t,cH,x)\ 
w{t,x)=\N{t,c\x)= , 

y\/{t,c^t,x)J 

where we separate the high oscillations from the slow time dependency of the 
solution. Here, W(t, 6, x) is a function defined on M x T x T*^ and taking values 
in C^. The equation for W is the transport equation (see (2.7)) 

(3.1) idt\N + ic^deVJ = -c(V)cW + c(V)-^F(W), 
with boundary condition 

(3.2) W(0,0,x) = •0(x). 

In a first step, we expand the previous equations with respect to the small 
parameter in the following sense: For a given k, we can write 



n>2 



^ n>2 

(3.3) = c2 + l|A;p + ^a„+ic-2"|fc'2"+2 

n>l 



AP SCHEMES FOR THE KG EQUATION IN THE NON-RELATIVISTIC REGIME 7 



for some coefficients a„ defining a power series of convergent radius 1. How- 
ever, as k is an arbitrary integer, the meaning of the previous expansion in 
terms of operators has to be understood in the sense of asymptotic expansions; 
this means that for a given sufficiently smooth function W, we can write 

c(V)cW = c^W - ^ AW + ^ a„+ic-2'^(-A)'^+iW, 

n>l 

in the sense that for ah € N, there exists a constant Cn such that if 
w G i7^+2A^+4 X i/«+2Af+4 for some s G N\{0}, then for all c> 0, 



1 ^ 

\\c{V),w - c^w + -AW-Y, an+ic-2"A"+i«;||^, < Cnc-^''-^\ 

n>l 

This last equation can be rigorously proved by using Taylor expansions of the 
function x i— t- t/1 + x > and Fourier decomposition of W. Similarly, we 
can write 

c(v)-i = (i - ^y'^" = 1 + + J;/3„c-2"(-A)^ 

n>2 

for some coefficients /3„. Hence in view of (2.6) the initial condition ip can be 
expanded in the sense of asymptotic expansion as 

^ = ^c-2>„(x), 

n>0 

with 

(3.4) ^0 ^i = -^Ar_], and ^„ = /3„(-A)" T ) . 

Similarly, Equation (3.1) can be expanded in the sense of asymptotic ex- 
pansion for smooth W as 

c'^iide + 1)W = - idtVJ + ^AW + F(W) + ^ c-2'^F„(W), 

n>l 

where 

F„(W) = a„+i(-A)"+iW + /3„(-A)"F(W). 
For later use we will give the first term Fi{\N): For W = (U, V)-^, 

AA2U + iA/(i(U + V))\ 



\A2v + iA/(i(V + U)), 



(3.5) Fi(W) = 

Note that for s > d/2, can prove the following estimates: 
Vn>l, ||F„(W)||^, <C„||W||;^ 
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where r is the polynomial degree of / (see (2.9)), and for some constant C„ 
depending on n and s. We now look for an asymptotic expansion for W as 

W(t,0,x) = ^c-2"W„(t,e,:E). 

n>0 

With this expression, we can expand the nonlinear term F(W) using the fact 
that / is polynomial. We obtain an expansion of the form 

F(W) = F(Wo) + c-'"dF(Wo) • W„ + 5„(W, | j < n - 1), 

n>l 

where S'n(Wj | i < n — 1) denotes a nonlinear term depending only on Wj for 
j < n — 1. Also note that this term is polynomial in the functions Wj. Here, 
di^(Wo) is the differential of F at Wq. It is a polynomial term of oder r — 1, 
in Wq if / is polynomial of degree r. 

Hence we get a recurrence relation of the form 

(3.6) (i5e + l)Wo = 0, and (i5e + l)Wi = (-i^j + iA)Wo + F(Wo), 

and for n > 2, 

(3.7) 

{ide + 1)W„ = i-idt + ^A)W„_i + dF(Wo) • W„_i + i?„(W,- 1 j < n - 2), 

where i?„,(Wj \j < n — 2) denotes a term depending only on Wj for j < 
n — 2 that is polynomial in the functions Wj and their derivatives. Using 
the previous estimates we obtain that for all n and all collections Wj, j = 
0, . . . , n — 2 of sufficiently smooth functions, 

(3.8) \\Rn{^j \j<n-2)\\^^< C„( _max ||Wj||^,+,,„_^, )^ 

where the constant C„ does not depend on the Wj's (but depends on n and 
s). Eventually, the initial condition (3.2) yields the conditions 

(3.9) Vn>0, W„(0,0,x) = V'„(x), 

where xjj^ is given by the expansion (3.4). 

To solve the previous induction relation (3.7) and (3.9), we impose the initial 
condition Wj = for j < 0. We also introduce the following space: we say 
that G{t,0,x) G %^ if G is a trigonometric polynomial in 9 with coefficients 
in C^{{0,T), X H^). In other words, we can write G as 

Git,9,x)= e*"V'^)(i,x), 

aSZ, |a|<p 

where p < +oo and with g^"'\t,x) G C^{{0,T), H'^ x H^). By construction, 
as the nonlinearity is assumed to be polynomial, it is clear that if for j = 
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0, . . . , n — 1, \Nj, is an element of ■'\ then we have (see (3.8)), 

Rn{\Nj\j <n-2)e'H'T. 
We will use repeatedly the following result: 

Lemma 3.1 (Solution of the homological equation) 

Let T G M and s > 0. Assume that G{t,6,x) G is given. Then there 
exists a unique solution W*(t,0,x) € Ti^ satisfying the equations 

(3.10) iid0 + l)\N^{t,e,x) = G{t,e,x)- — [ e~''^G{t,e,x)de, 

2vr J J 

and 

(3.11) VtG(0,r), VxeT*^, — [ e-'^\N^{t,e,x)d9 = 0. 

Moreover, for all function w{t,x) £ C^{{0,T), H'^ x H^), the function 

\N{t, 9, x) = e'^w{t, x) + W,(i, x). 
satisfies the equation (3.10). 

Proof. — We decompose G(t, 6, x) in Fourier in 9: we can express this function 

as 

G(t,e,x)= ^e-V'^)(i,x), 

|a|<P 

with p < +00 and where g'^''\t, x) £ C^((0, T), H" x H""). Note that 

^ I e-''G{t,9,x)d9 = g('\t,x). 

Let us seek W^, under the form 

\N4t,9,x) = J2 e'"^ti;l"^(t,x), 

for some unknown functions wi'^^ in C^{[0,T), x H^). Equation (3.10) can 
be written as 

V|a| <p, {l-a)wi''\t,x)=g(''\t,x)-6lg^'\t,x), 

where (5^ denotes the Kronecker symbol. It is clear that this equation is solv- 
able for the functions wi^\t,x), a ^ 1 and that with the condition (3.11) we 
can set {t, x) = 0. The last statement is a consequence of the fact that 
the kernel of the operator id^ + 1 is the space spanned by e*^. □ 

We are now ready to state the main result of this section: 
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Theorem 3.2. — LetNen,s> d/2, cq > and assume that ip G i7s+27V+4 
and 7 G f{s+2N+A ^j,g giy^^i. Then there exists T = T{s, N) > 0, and a se- 
quence of functions \Nn{t,6,x) G 'H^''^^^ ^\ satisfying the following proper- 
ties: Wo = wqc^^ where wq satisfies the equation 

(3.12) idtWo = lAwo + {F){wo), too(O) = V'o> 
where ipQ is defined in (3.4) and where by definition, 

(F)(w) = — / I a,/, : ] d6, if w = [ 

^ ^ 2TrJo \f{v + e-^'^u)J \v 

and for all n, \Nn solve the equations (3.7)-(3.9) for n = 0, . . . , N , and can be 
decomposed into trigonometric polynomials of the form 

(3.13) W„= E 

l«l<Pn 

with pn < +00, and where wl?'\t,x) G C((0, T), x /7«+2(n-Af)). 
Moreover, there exists a constant Cn such that for c > cq, if w{t) denotes the 
solution given by Proposition 2.2 with initial value ip given by (2.6) in terms 
of if and 7, then 

N 

(3.14) Vt<r, \\w{t,x)-Y,c~^"^n{t,c-\x)\\^^<CNC-^^-^. 

n=0 

Proof. — Let us consider the recursion (3.7) together with the initial condition 
(3.9) where the ■0^ are given by (3.4). For n = 0, we obtain the equation 

i{d0 + l)\No{t,9,x) = 0, and W(0, 0, x) = t/>o(x) 
which imphes, according to the previous lemma, that 

(3.15) \No{t,e,x) = e'^w'l^\t,x), and wl^^\o,x) = ipQ^x) 

that is we have determined = for a 7^ 1 (and set po = !)• 
For n = 1, the equation (3.6) is written 

(3.16) i{de + l)Wi(t, d, x) = i-idt + ^A)Wo + F(Wo). 

According to the result of given in the previous lemma (see (3.10)), this equa- 
tion will be solvable if the right-hand side satisfies 

/ i-idt + i A)e-^^Wo(t, e, x) + e-^^F(Wo(t, 9, x))d0 = 0. 

In view of (3.15), this equation can be written as 

{idt-'^A)w^^\t,x) = ^ [ e-''F{e''w^^\t,x))de, 
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which gives the equation (3.12) with the boundary condition (3.9). Now let T 
be such that WQ{t, x) is a solution to the previous equation with initial value 
w{0,x) = ipoix) in the space C\{0,T), H'+'^^+^ x H'+^^+^). Then we can 

set w'^\t,x) = WQ{t,x) and solve for the functions w^^\ a / 1 according to 
the previous lemma. Note that pi is given by the degree of the polynomial in 
the nonlinear right-hand side. 

Now assume that n > 1 is given, and that w^j^\t,6, x), a < pk are known 
functions in C^((0, T, _H's+2(Af-fc)+4 ^ ^s+2{N~k)+i^ for A; < n - 2 and pk < oo, 

and that w^}^-^^{t,9, x) is determined for a ^ 1. 

Consider Equation (3.7). Using (3.8) and the induction hypothesis, 
the right-hand side of this equation belongs to C^{{0,T, H'^^'^^^~"'^'^'^ x 
^,+2(7V-n)+2)^ and is a polynomial in e*^ of degree pn > Pn-i, because the 
nonlinearity is a polynomial. To be solvable, the right-hand side of (3.7) has 
to be orthogonal to e~*^ in L^(T). By induction hypothesis, only the term 

w^}_i has to determined. In view of (3.7), we thus see that w^}_i has to satisfy 
a linear equation of the form 

(3.17) {-idt + lA)w^^l^{t, x) = m„_i(t, x) • wj^^ + r„_i(t, x), 

where mn-i{t, x) is a linear multiplicative operator with coefficients in the 
space C^{{0,T),H'+^^+^), and where r{t,x) G Ci((0, T, /7^+2(^-")+2 x 
jjs+2{N~n)+2y "Yhe initial condition associated with this equation is written, 
according to (3.9), as 

wl'l,{0,x) = 'il,^{x)+ ^2i(0'^)' 

H<Pn-l, a^l 

where tpni^) is given by (3.4) in terms of <p and 7, and where the other terms 
in the right-hand side are known by induction hypothesis. 

As the equation (3.17) in linear with coefficients bounded on the time inter- 
val (0, T); we deduce that w^^l^{t, x) is weh defined in C^((0, T, H^+2{N-n)+2 ^ 
jjs+2{N-n)+2y ^sing the previous Lemma, we can then solve for the terms 
wi^^ for a / 1 in Ci((0, T, /i'^+2(^-")+2 x i/^+2(W-n)+2). 

To prove the asymptotic expansion, we note that by construction the term 

AT 

w^{t,x) := ^c-2«W„(t,c-2t,a;) 

n=0 
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belong to C^((0, T, H^^'^ x H^^'^) and satisfies the equation (2.7) up to an error 
or order c"^^^^^-*, that is 

(3.18) VtG(0,r), 

idtw^{t,x) = -c{V)cW^{t,x)+c{V)-^F{w^{t,x)) + c-^^-^r^{t,x), 

w^{0,x) = -0(2^) + c-^^-^Vf 
with 

r^{t,x) £C\{0,T,H' X H'), and tf^^ (x) £ H' x H' , 
satisfying the following uniform bounds: for a given cq > 0, there exists a 
constant Cat such that 

Vc > CO, Vt G (0,r), \\w''it,x)\\^^ + ||rf (t,x)||^, + (x)||^, < C^v- 

Using Proposition 2.2, we can assume that the exact solution w{t) is also 
uniformly bounded (with respect to c) on (0,T) in H^. We conclude by using 
(3.18) and standard comparison arguments based on the Gronwall Lemma on 
(0, T) and in with s > d/2. □ 



4. First terms 

The previous result implies the existence of an asymptotic expansion of 
u{t,x) and v{t,x) solutions of (2.4)-(2.5) up to any arbitrary order, provided 
the initial condition is smooth. Note that we will use the following notations: 
for all n > 0, W„ = (U„, V„)^ with, according to the decomposition (3.13) 

(4.1) U„= ^ e'^%i^Ht,x), and V„ = ^ e''^%i-\t,x). 

|a|<p„ |a|<Pn 

This implies an asymptotic expansion of the solution u{t), v{t) and z{t) of the 
form 

u{t, x) = C~'^"'Un{t, C^t, x), v{t, x) = C^'^"'Vn{t, C^t, x), 
n>0 n>0 

and 

z{t,x) = ^c"2"z„(t,c2t,x), 

n>0 

with the relation Zn = ^{un + v^)- Note that here, the equality signs must be 
understood in the sense of asymptotic expansion. Remark also that in the case 
where (p and 7 are real, that is when the solution z{t) is real, then we have for 
all n, Un = Vn- Using the result of Theorem 3.2, ZQ{t,x) is an approximation 
of the solution z{t, x) over the time interval (0, T) with an error order 0{c~'^), 
while 

(4.2) zo{t,x) + i,zi{t,x) 
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approximates the exact solution z{t, x) over the time interval (0, T) up to terms 
of order c~^. The goal of this section is to express the two first terms zq and 
zi . For simplicity we consider the one dimensional case and assume the initial 
values ip, 7 to be sufficiently smooth so that Theorem 3.2 can be applied. Note 
that a generalization to higher dimensions follows the line of argumentation. 
We will use these computations to derive our numerical schemes in Section 6. 
For reasons of ease we will focus on the linear case f{z) = Xz, see Section 4.1, 
and a cubic-nonlinearity f{z) = X\z\'^z, see Section 4.2. Nevertheless we state 
the equations for Ui and Vi in the general case f{z) = X\z\'^Pz, p £ N. 

The first term zq approximating the exact solution of Klein-Gordon equation 
(1.1) up to terms of order is easily derived from Theorem 3.2 using (2.3) 
and (2.6). More precisely, we have 

zo{t,x) = ^e**'='no(t,x) + le-'^^^v^it, x), 

where uo,vq solve the nonlinear Schrodinger system (3.12). 

In order to explicit the term zi{t, x), we expand the non-linearity /(^(U+V)) 
in a Taylor-series, up to terms of order c^^. With the notation Z^. = ^(Ufc+V^) 
and keeping in mind that f{z) = f{z,z) = X\z\'^^z = Xz'^^^z^, we obtain 

/(i(U + V)) = /(Zo + ^Zi) + 0(c-4) 

= /(Zo) + ^ {dJ{Zo, To)Z, + aj/(Zo, Z^)ZT) + 0(c-4) 

= /(Zo) + ^[{P+ l)|Zo|2?'Zi +p|Zo|2f-2(Zo)2Zl] + 0(c-4). 

Furthermore, in the above notation we calculate that 

A/ = d,j{d,z)^ + 2d,^f{d,z){d,z) + dj{d^^z) + ckzf{dj.f + d^f{d,,z), 

where in this formula / is evaluated in (Z, Z). 

Thus, we obtain by (3.7) with n = 2 (using (3.5) and the above calculations) 
that the homological equation in terms of U2 is written 
(4.3) 

{ide + 1)U2 = - idtUi + ^AUi + IA^Uq 

+ A(p + l)|i(Uo+V^)pfi(Ui + VT) 

+ Ap|i(Uo + v^)P^-2(^(Uo + V^))2^(Vi + U7) 

+ Xip + l)|i(Uo +\^)pPiA(Uo + V^) 

+ Ap|i(Uo + Vo)|2f-2(l(Uo + Vo))2|A(Vo + Do) 

+ Aip(p+l)|i(Uo + V^)|2p~2i(u^ + Vo)(94(Uo + V^))2 

+ xp{p + i)|i(Uo + v^)|2p-2i(Uo + v^)|a4(Uo + y~o)\^ 

+ Aip(p-l)|i(Uo + V^)|2p-4(|(Uo + V^))3(a4(L^ + Vo)), 
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and a similar equation for V2 with Ufc replaced by V^, k = 0, 1,2, and vice 
versa. 

4.1. Formula for the first terms in the linear case. — In the linear 
case we have f{z) = \z. In this case we obtain the following limit equation 
for uq and vq 

^^^^ {-idt + lA + ^)uo = 0, no(0) = (/J - i7, 

{-idt + lA + ^)vo = 0, vo{0) = ^ + ij, 

see (3.12). For Ui we obtain by (3.16) 

{ide + l)Ui = i-idt + lA)uoe'<> + |A(uoe*^ + v^e-'^). 

Thus, plugging the Fourier expansion (4.1) of Ui in 9, we obtain with the 
orthogonality relation (4.4) for uq that 

Thus, n^"^ = for a / ±1 and u[ = \Xvo, and we set ^1 := u^i \ so that we 
can write 

(+ a ^\ — ^ 



(4.5) Ui(t, 6, x) = IXv^it, x)e-^^ + ^t, x)e'' , 



Note that similarly \/i{t,6,x) = \Xuo{t,x)e~^^ + r]i{t,x)e^^ with rji := v\^\ In 
order to determine .^i and rji we analyze (4.3) in the linear case, i.e. p = 0, 
which reduces to 

{ide + 1)U2 = -idt(Ji + iAUi + gA^Uo + A^(Ui + V^) + |AA(Uo + 

To be solvable, the right-hand side of the equation needs to be orthogonal to 
e~'^ which yields the following equation for ^i{t,x): 

i-idt + 5A + ^)ei = -gA^no - lA^uo - IXAuQ. 

Using twice the Schrodinger equation (4.4) for uq, the above equation can be 
rewritten as 

{-idt + 5A + |)ni = IdttUf). 

The same equation holds true for rji with ^1 and uq replaced by r]i and vq, 
respectively. The initial values ^i(O) and r]i{0) have to be chosen such that 
Wi(0,0,x) = tpiix), see (3.9). As Ui(0) = \Xv^{0) +6(0) and Vi(0) = 
|Auo(0) + r/i(0) we thus obtain by (3.4) the following initial value problems 
for 6 and rji 

{-idt + lA + I) Ci = ^duuo, ^i{0) = liAuo{0)-{A + X)v^{0)), 
{-idt + lA + ^)rji = ldttVo, r?i(0) = i(At;o(0)-(A + A)n^(0)). 
We summarize our findings in the following corollary. 
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Corollary 4-1 (First terms in the linear case). — In the linear case, 
i.e. f{z) = Xz, the first and second terms in the asymptotic expansion (4.2) 
read 



(4.7) zi{t,x) = |(no(t,x)e^^'*+t)^(t,x)e-^^'*) + i(a(t,x)e*^'*+r?r(t,x)e-^^'*), 



where uq, vq and^i, rji solve (4.4) and (4.6), respectively. Then the functions 
ZQ{t,x) and ZQ{t,x) + ^zi{t,x) approximate the exact solution of (1.1) in 
up to terms of order and c~^ for initial values in H^'^'^ and H'^^^ , 
respectively. 

Note that both corrections terms (4.4) and (4.6) can be solved exactly in 
Fourier. Thus, in the linear case we can state the correction terms zq and zi 
in an exact way. See Section 6.1 for further details. 

4.2. Formula for the first terms in the cubic case. — We consider 
now the case of cubic nonlinearity, i.e. p = 1 with the previous notations. 
Moreover, for simplicity we consider the case, where uq = vq, i.e. z{0) and 
dtz{0) £ M. Thus, /(i(Uo + Vo)) = ^{uoe'^ +u^e-'^f and the limit equation 
for Uq reduces to 

(4.8) - iSfUo + jAuo + |A|nopno = 0, tto(O) = v? - 27, 

see (3.12). Plugging the Fourier expansion (4.1) of Ui in 9 into (3.16), we 
obtain with the orthogonality relation (4.8) for uq that 



Hence, u\ = for a ^ —3, —1, 1, 3 and 

u^f^'^ = ^{u^f, u[-^^ = %\uo\^u^, and uf^ = -^{uof. 



zo{t,x) = le''^''^uo{t,x) + '^''^vo{t,x) 



and 




Setting ^1 



u 



(1) 



we obtain 



1 ' 



(4.9) Ui 



-3ie 



+ eie^'(= Vi). 
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In order to determine the function we need to state the equation for U2 = 
X^aez ^2'*^^*"^' which reduces in our cubic setting to 

+ ^a((uo)V^^ + 2\uo\^ + (n^)2e-2^^)(A^xoe^' + An^e"^^) 

see (4.3). To be solvable, the right-hand side of the equation needs to be 
orthogonal to e~*^ which yields the following equation for ^i{t, x) = u^i\t, x): 

{idt - iA)ei - A||uoP6 - Xliuo)Ti = lA^no + \^§-,\uo\\o + lAfi^uo). 

The initial value ^i(O) needs to be chosen such that Ui(0) = — IA7, see (3.9) 
as well as the definition (3.4). Thus, by (4.9) we obtain that 

6(0) = jsMO? - ^miOf - iko(0)|2u^(0) + |A(no(0) - mio)). 
Using (4.8) we rewrite the equation for ui as 

(ia^ - iA)6 - AfhoPa - A|(uo)'6 

(4.10) = lidtAuo + A^^lnol^no + A^A(|moPmo), 
6(0) = ^^0(0)=^ - i^MOf - iko (0)1 2^2^(0) + ^A(tio(0) - MO))- 

Again we summarize our findings in the following corollary. 

Corollary 4-2 (First terms in the cubic case). — For a cubic non- 
linearity, i.e. f{z) = \\z\'^z, and for real initial values functions 2(0,2;) = 
(/^(x) G M and dtz{0, x) = (?^{x) G M, the first and second terms in the 
expansion (4.2) read 

zo{t,x) = 2e**''^no(t, x) + 4e~**''^n^^(i, x), 

and 

zi{t, x) =§|no(t, x)\^{uo{t, x)e*'='* + uoit, x)e-*'='*) 

(4.11) - A (^uoit, xfe^''^'' + n^(t, xfe-^^'^'') 

+ i(6(t,x)e-'* + 6(t,x)e--'*), 

where uq and^i solve (4.8) and (4.10), respectively. Then the functions ZQ{t, x) 
and zo{t,x) + ■^zi{t,x) approximate the exact solution of (1.1) in up to 
terms of order and for initial values in H^'^^ and H^~^^, respectively. 
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5. Energy and charge density conservation 

Analyzing the conserved quantities of the Klein-Gordon equation (1.1), we 
follow Dirac's approach via the "hole" theory, i.e. antiparticles are related 
to negative energy eigenstates (n particle with positive charge, v antiparticle 
with negative charge). This allows us to overcome the problem of indefinite 
probability density (i.e. if we would consider z{t) as a probability amplitude, 
its probability density p{z) = {zdtz — zdtz) which is naturally defined by 
the continuity equation, would allow positive and negative values.) For further 
details on the quantization of the Klein-Gordon field we refer to [4]. 

In this interpretation {p{z) being the c/iargfe-density) the conserved physical 
quantities of the Klein-Gordon equation are the charge 

(5.1) Q{z) = ^ j zdtz - zdtzdx = J Re ^-J-dtzz^ dx, 
as well as the energy 

(5.2) E{z)= [ \ldtz\^ + \Vz\'^ + \cz\^--^\z\^P+^dx. 

Jt 

Hence for all times t where the solution exists and is smooth enough, 
Q{z{t)) = Q{z{0)) = [ Im(^7)dx, 

and 



T 



E{z{t)) = E{zm = cmMl, + hlliO + \V^\' - ^M''"-' drr. 

Using the result given by Theorem 3.2, and under smoothness assumptions 
on (fi and 7, we can easily prove that the charge and energy admit asymptotic 
expansion with respect to c~^. In this section, we give the expression of the 
first terms of these expansions. We commence with the asymptotic expansion 
of the charge Q{z), i.e. 

Q{z) = Qo(Uo, Vo) + ^c-2'=Qfc(Uo, . . . , Ufc, Vo, . . . , V^). 

k>l 

Using (2.2) we can easily (as we are working on the torus T and thus integration 
by parts yields no boundary-terms) rewrite the charge as 



(5.3) 



Q(z) = ^j^Re(^{^^iu-v))iu + v)^ dx 

= \ I h|2-|z;|2 + -L(|Vn|2-|V^;|2)drc + 0(c-^) 
4 Jf Ac 
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Thus, using the asymptotic expansion given in Theorem 3.2, we obtain 

(5.4) go(Uo,Vo) = ^(||^xo|li2-|bo|li2), 

where uo,vo solve the nonhnear Schrodinger system (3.12) and hence in par- 
ticular their norms are conserved. Thus, the first term of the charge reads 

(5.5) Qo(Uo(t),Vo(t)) = \ (llno(O)lli. - |bo(0)||i.) . 

Using the relation z = ^{u + v) as well as dtz = ^ic{V)c{u — v), see (2.4), we 
obtain 

E{z) = l [ c2(|n|2 + |^;|2) + |Vn|2 + |V^;|2--^|i(n + i;)|2p+2dx + 0(c-2). 

Hence similarly to Q{z), we obtain an asymptotic expansion of the form 
Eiz) = c^E_2i^o, Vo) + Yl c"''^fe(Uo, . . . , Ufc, Vo, . . . , Vfc), 

fc>0 

with 

^_2(Uo,Vo) = ^(||uo||i2 + ||?;o|li2). 

In order to analyze the asymptotic behavior of the energy E{z), we first 
have to subtract the relativistic rest-energy c^Qmass (see below) as from the 
energetic point of view the non-relativistic limit regime means that £ := E — 
c^Qmass is Small. According to (5.3), we can write 

Q{z) = Q^{u) + Q^v), 

where in the sense of asymptotic expansion (see (3.3)), 

Q"(n) = ^ / + JL|Vn|2 + ^aj-c-2j|V^u|dx, 

•^T j>2 

and where Q^{v) denotes the part of Q{z) which contains the antiparticle 
component, i.e. 

With this notation, we define the rest energy c^Qmassiz) = 2c^{Q'^{u) — Q^{v)). 
Using these notations, we obtain an asymptotic expansion of the form 

£{z) := E{z) - 2^{(r - QT) 
(^•^^ = fo(Uo, Vo) + c-2^ffc(Uo, . . . , Ufc, Vo, . . . , Vfc). 

k>\ 

Using the asymptotic expansion given in Theorem 3.2, we calculate that the 
first-order correction term of the energy if (2) in (5.6) is given by 

2p+2 



(5.7) 



^^o(Uo,Vo)= /^(|VnoP + |V7;oP)-^ 
Jt 4 ' ' p+1 



i(no+we-2-'*; 



dx. 
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where uq,vo solve the nonhnear Schrodinger system (3.12). 



6. Construction of robust numerical schemes 

The aim of this section is to show that we can use the derived asymptotic ex- 
pansions to construct numerical schemes for the non-relativistic Klein-Gordon 
equation (1.1) which are robust with respect to the speed of light c. Firstly we 
underline our theoretical results in the linear case, where the exact solution 
to (1.1) is analytically known, see Section 6.1. Finally, in Section 6.2 we con- 
struct robust numerical methods for the cubic non-relativistic Klein-Gordon 
equation and discuss their convergence. Numerical experiments are included. 



6.1. Linear case. — In the linear case, i.e. p = 0, our model problem reads 

\dttz - Az + c^z = Xz, 
c 

^^■^^ z{0,x) = ^ = Y,ipae'--, 9,z(0,x) = c27 = c2 5^W" 

with the exact solution 

z{t,x) = (^ipa cos{ct\/ a? + — A) 
(6.2) ""^^ 

+ , ^ \ ^ 7asin(ctA/a2 + c2-A))e^°^ 
V + — A ^ 

In the linear case we can solve the first- and second-order limit systems exactly 
in Fourier, see Corollary 4.1 and equation (4.4) and (4.6), respectively. Thus, 
we can state the first- and second correction terms zq and zi in an exact way. 
A simple computation shows that 
(6.3) 

zo{t,x) = Y, ((^acos ((c2 + ^(a2 - A))t) +7asin ((c2 + \{a^- A))t))e*'^^ 

Thus, we nicely see that the difference between the exact solution z of the 
linear Klein-Gordon equation and its first-order correction term zq in the non- 
relativistic limit c ^ 1 are the formal approximations 

CA/a2 + c2-A«c2 + -(a2-A), , ^ 1, 

2^ Va2 + c2 - A 

which are of order 0{c~'^) for sufficiently smooth initial values and 7. 
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For the second correction term zi we get after a short computation 

Zi{t,x) =^ \(pa "g^-* * sin{Ua,ct) 



(6.4) 



where 



+ 7a((-^^V^)sin(a;a,ci) - ^" g^-* * cos(a;a,ct) 



Wa,c = + \{a^ - A). 



Again this can easily be motivated by formally taking the next approximations 
in the exact solution z, i.e. 



+ c2 - A « c2 + i(a2 - A) - -^(a^ - A)^, 



and 



1 



2c2 ' 



Va2 + c2 - A 

which are of order 0(c~^) for sufficiently smooth initial values and 7. 

The numerical results for the initial values 

{2 + i) , , (1 + i) . , . 1 , . 

^= ^ cos(x), 7 = ^ sm(x) + - cos(x) 

and A = — 1 on the time interval [0, 1] are given in Figure 1. 



I 

-^10-' 
-ft 10-' 




10" 



10 



10' 



'2c 



Figure 1. Numerical simulations for the linear Klein-Gordon equa- 
tion: Difference between the exact solution z and the first- and 
second-order approximation terms y = Zq (red, cross) and y — 
zq + c~^zi (blue, star). The slope of the dash-dotted and dotted 
line is two and four, respectively. 
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6.2. Cubic non-linearity 

problem reads 

1 



In the cubic setting, i.e. p = 1, our model 



(6.5) 



duz-Az + c''z = X\z\''z, z{0) = ip, dtz{0) 



c^7 



and the limit Schrodinger system (1.3) reduces to 



(6.6) 



idtuo 
idtvo 



Auo 
-AvQ - 



A 



A 



|no|^ + 2|wo|^)no 



l^ol^ + 2|no|^)uo 



= 0, 
0, 



no(0) 
vo{0) -- 



If - Z7, 



((5 + Z7, 



-ic't 



i.e. the first-order correction term zq is given by 

(6.7) zo{t) = iuo(t)e*^'* + lv^{t)e 

where (uq, fo) solves the system of nonlinear Schrodinger equations (6.6). Here 
the enhancement of our Ansatz is clearly seen: the limit Schrodinger system 
can easily be simulated without any c-dependent time step restriction, and a 
well suited approximation to the exact solution in the non-relativistic limit 
regime is then obtained by simply multiplying the numerical approximations 
of uq and vo with the highly-oscillatory phases e^^ * and e~*'^ *, respectively. 

For the numerical time integration of the first-order limit system, i.e. the 
Schrodinger system (6.6), we apply an exponential Strang splitting method, 
where we split the limit system in a natural way into the kinetic 



(6.8) 

and potential part 
(6.9) idt 



idt 



Uq 


= 1a 


Uo 


Vo 


2 


Vo 



Uo 


A 


"koP + 2|woP 







Uq 


.^0. 


~ 8 





|woP + 2|'UoP. 




Vo 



The exact flow $*(no(0), t;o(0)) = ^ip_^_p{uo{0),vo{0)) of the coupled system 
(6.6) at time = nr is then approximated by 



(6.10) w (^'p'^o^^o^p'- 

where $^(ito(0),wo(0)) and ^'*p(uo(0), vo(0)) are the exact flows associated to 
the kinetic Hamiltonian (6.8) and potential Hamiltonian (6.9), respectively. 

Note that the kinetic equation (6.8) can be solved exactly in Fourier space. 
Furthermore note that in the potential equation (6.9) the moduli of uo and vo 
are conserved quantities, i.e. |mo(*)P = l^o(0)P an bo(*)P = l^o(0)P for ah t 
and thus we can solve both split equations in an exact way. 

For the space discretization, we use pseudo-spectral Fourier with mesh-size 
h and grid points Xj = 2j/i7r, j = —K, . . . ,K — 1 where K = 1/h £ N is the 
number of frequencies. Note the fully discrete scheme can be easily implement 
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using the Fast Fourier transform algorithm. We refer to [17, 9] for convergence 
results over finite time of such splitting schemes. 

Example 6.1 (First-order approximation term). — In the first exam- 
ple we are interested in the asymptotic expansion up to the first-order cor- 
rection term zq. The asymptotic expansion derived in the previous section 
allows us the following convergence result for the natural choice of numerical 
approximation to zq. 

Theorem 6.2. — Let s > d/2 and assume that ip andj are smooth functions 
defined on the one- dimensional torus T. Then there exist constants T, cq, C, Hq 
and To such that the following hold: for c > cq, t < tq and h < h^, if we define 
by 

(6.11) z;;'^ := l^'^^^'*" + i^^e-*^'*" 

the numerical approximation of the first-order approximation term ZQ(t) at 
time tn = riT, where Uq''' and v^'^ denote the numerical approximations to the 
first-order limit system (6.6) obtained by the Strang splitting (6.10) at time 
tn = nr with a discrete grid of mesh-size h; then if z{t) denotes the exact 
solution of the equation (6.5), we have 

\\z{tn) - ^o'^IIl^ < C{t^ + h' + c-2) for all t„ < T. 

Proof. — We have 

(6.12) \\z{tn) - Zq'^\\l2 < \\z{tn) - Zo{tn)\\L2 + \\zo{tn) " Zq'^\\l2. 

By Theorem 3.2 we have over the time interval (0, T) and c > cq, 

(6.13) Mtn) - Zo{tn)\\L2 < Cc-\ 

To bound the second term we use the following inequality 

< ||'"o(*n) - Uq'^Wl^ + \\vo{tn) " 'yo'''||L2- 

Using arguments similar to [9, Theorem IV. 17] (see also [17]), it is then pos- 
sible to prove that 

(6.14) ||no(t„) - u^'''\\l2 + \\vo{tn) - t^o'^lli^ < C{t^ + h'), 

if the exact solution is smooth enough (i.e. belongs to some Sobolev space H^' 
with s' > s sufficiently large) . Note that the proof of the fully discrete scheme 
in [9] is given only for the Lie-splitting, real initial data (that is (p and 7 real 
valued functions) and uses function spaces based on the Wiener algebra (the 
il spaces). However, the rigorous proof of (6.14) is essentially a variation on 
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the same theme and does not present any particular difficulty as soon as the 
solution is smooth enough. Using this bound, we obtain that 



(6.15) \\zo{tn)-Z^'''\\L2<C{T^ + h'). 

Plugging the bounds (6.13) and (6.15) into (6.12) yields the desired conver- 
gence result. □ 

The numerical results for the initial data 

(2 + i) , . {l + i) . . 1 . ^ 
ip = ^ cos(x), 7 = ^ sm(x) + - cos(2;) 

and A = — 1 on the time interval [0, 0.1] are given in Figure 2. As a reference 
solution for z{t) we take the numerical approximation of the Klein-Gordon 
equation (6.5) obtained with an exponential Gautschi method, proposed in 
[2], with a very small time step r satisfying the necessary condition r < (?h. 
Note that the convergence result for the exponential Gautschi method given 
in [2, Theorem 9] shows a convergence rate of order 0(r^c^ + h^) under some 
a priori smoothness assumptions and uniform bounds for the exact solution. 



5b 
I 



ft 



10" 




10 10 

\/2c 



Figure 2. Numerical simulations of Example 6.1: Difference be- 
tween the reference solution z obtained with an exponential Gautschi 
method with time step r = 10~^ and the first-order approximation 
solution y = Zq''' obtained with r = 10~^. The slope of the dash- 
dotted line is two. 



Example 6.3 (Second-order approximation term) 

In this example we simulate the first- and second-order approximation terms 
zq and zi given in (6.7) and (4.11), respectively. For the numerical approx- 
imation of the first approximation term we proceed as in Example 6.1. For 
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the numerical approximation of the correction term zi, we solve the equation 
for ^1, i.e. equation (4.10), with a symmetric Strang splitting based on the 
following split-equations: The kinetic equation 

which can be solved exactly in Fourier, and the linear potential equation 
(6.16) idt^ = AfluoPC + A|(no)2e + |A2no + X^§^\uo\^uo + X^A{\uo\uo), 



= Im(^^o), 

A^^|uo|^tto + XjqA{\uo\uo). Then we can rewrite the potential equation as 
(6.17) 



which we solve as described in the following. We set ao = Re(no), / 
Q = Re(^), /3 = Im(^), and furthermore go{t) = g{ao{t), (3Q{t)) 



dt 
We set 



2ao/3o 
-((/3o)2 + 3(ao)') 



((ao)' + 3(/3o)') 
-2ao/3o 





a 


+ 


Im(c/o) 




P. 


-Re(5ro) 



A(ao,/3o) = -A^ 



2ao/3o 
-((/3o)2 + 3(ao)' 



-2ao/3o 



and approximate the potential equation (6.17) with an exponential trapezoidal 
rule, i.e. for afc(t„), ^ Mtn), /c = 0, 1, = g{a^,f3^), t„ = nr, we 

set 



(6.18) 



,i(AK+i,/3o"+^)+AK,/3y)) 



a" 



+ 



Im(ffo") 
-Re(<7o") 



+ 



Im«+^) 



The numerical approximation ^"+^ to the exact solution ^(t„+i) of the po- 
tential equation (6.16) starting in a given point is then simply given by 



tra+l 



a 



n+l 



n+l 



Concerning the space discretization, we use again the Fourier pseudo- 
spectral method described in the previous Section, and we denote by h the 
corresponding mesh size. Note that in this setting, the approximation of 
can be easily made using FFT. 

Let Uq'^ be defined as in the previous section, and let ^"''^ denote the nu- 
merical approximation to the full system (4.10) at time t„ obtained by the 
fully discrete Strang splitting described above, where we make half a step 
with the kinetic, a full step with the potential and again half a step with the 
kinetic part. Then, as the numerical approximation to the exact second-order 
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approximation term at time t„ = we take 

(6.19) Z^''^ + C-'z'^''' := i(l + -^3A|^n,h|2)(^„^.c^*„ ^^^-.c^t„ 



For this approximation we have the following convergence result. 

Theorem 6.4- — Let s > d/2 and assume that Lp andj are smooth functions 
defined on the one-dimensional torus T. Then there exist constants T, cq, C, Hq 
and To such that the following hold: for c > cq, t < tq and h < Hq, then if 
z{t) denotes the exact solution of the equation (6.5), we have 

Mtn) - izo'"" + c-'2"''^)||l2 < C7(r2 + h' + c-^) for all tn < T, 

where the approximation z^'^ + c^^z"'^ is defined in (6.19). 

Proof. — Note that by Theorem 3.2 we have over the time interval (0, T) and 

c > Co, 

Mtn) - (Zoitn) + C-2zi(t„))||^2 < Cc'^. 

The rest of the proof follows the line of argumentation given for the proof to 
Theorem 6.2, since the scheme described in (6.18) is second-order convergent. 
We do not give the detail of this proof. Again it is a variation on the same 
theme as in [9, 17]. □ 

The numerical results for the initial data 

(p = cos(x), 7 = ^ sin(x) + - cos(x) 

and A = —1 on the time interval [0,0.1] are given in Figure 3. 

Example 6.5 (Energy and charge-density deviation) 

Here we simulate the deviation of the the first energy and charge correction 
term (5.7) and (5.5), respectively. Also, we simulate the energy and charge 
deviation of the first-correction term ZQ{t), i.e. E{zo{t)) and Q{zo{t)), where 
the energy and charge are defined in (5.2) and (5.1), respectively. 

We set E{zq'^) = E{uQ'^e'^'^^^" jV^'^e^'^^*") the numerical energy of the first 
correction term zq, for t„ = nr, and where Uq'^ and Vq''* are the numerical 
approximations of the first limit system (6.6) again obtained by the Strang 
splitting (6.10). Similarly, we set Sq''^ = £o{uQ'^e^'^^^" , v^'^e^'^^^") the approx- 
imation of the first term in the asymptotic expansion of the energy. In the 
same way, we define the numerical charge Q{zq'^) = Q(uQ''^e*'^^*", Ug'^^e*'^^*") 
of the first correction term, as well as the numerical approximation of the first 
term in the asymptotic expansion Qq'^ = Qo{uq'^ ,Vq'^). 

We carry out the simulations on the time interval [0, 1] with the time step 
T = 10"^. As the initial values we choose the same as in Example 6.1, but 
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10° 10' 10^ 



a/2c 

Figure 3. Numerical simulations of Example 6.3: Difference be- 
tween the reference solution z obtained with an exponential Gautschi 
method with time step r = 10~^ and the first- and second-order ap- 
proximation solutions y = Zq (red , cross) and y = Zq'^ + c"^^"''' 
(blue, star) obtained with r = fO~"^. The slope of the dash-dotted 
and dotted line is two and four, respectively. 

normalize them with respect to the norm. We nicely see that the first en- 
ergy correction term £q and the first charge correction term Qq are conserved 
by the Strang splitting, see Figure 4 and Figure 5, respectively. Further- 
more, Figure 6 and Figure 7 underline that \Q{zq''^) — Q{zq'^)\ = 0{c^^) and 
\E{z','')-E{z^'')\ = 0{l). 
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Figure 7. Numerical Simulation of Example 6.5: Energy deviation 



AP SCHEMES FOR THE KG EQUATION IN THE NON-RELATIVISTIC REGIME 29 



[9] E. Faou, Geometric numerical integration and Schrddinger equations, European 
Math. Soc (2012). 

[10] E. Faou, B. Grebert, Hamiltonian interpolation of splitting approximations for 

nonlinear PDEs, Found. Comput. Math. 11 (2011) 381-415 
[11] E. Faou, B. Grcbcrt, E. Paturel, Birkhoff normal form for splitting methods 

applied to semi linear Hamiltonian PDEs. Part I: Finite dimensional discretization, 

Numer. Math. 114 (2010) 429-458 

[12] E. Faou, B. Grebert, E. Paturel, Birkhoff normal form for splitting methods 
applied to semi linear Hamiltonian PDEs. Part U: Abstract splitting, Numer. Math. 
114 (2010) 459-490 

[13] J. Ginibre, G. Velo, The Global Cauchy Problem for the Non Linear Klein- 
Gordon Equation, Math. Z. 189 (1985), 487-505. 

[14] E. Hairer, C. Lubich, G. Wanner, Geometric Numerical Integration. Structure- 
Preserving Algorithms for Ordinary Differential Equations. Second Edition. Springer 
2006 

[15] M. Hochbruck, C. Lubich, A Gautschi-type method for oscillatory second-order 

differential equations, Numer. Math. 83 (1999), 403-426. 
[16] S. Jimenez, L. Vazquez, Analysis of four numerical schemes for a nonlinear 

Klein-Gordon equation, Appl. Math. Comput. 35 (1990), 61-94. 
[17] C. Lubich, On splitting methods for Schrddinger-Poisson and cubic nonlinear 

Schrddinger equations, Math. Comp. 77 (2008), 2141-)2153. 
[18] S. Machihara, N. Masmoudi, K. Nakanishi, Nonrelativistic limit in the energy 

space for nonlinear Klein- Gordon equations. Math. Ann. 322 (2002), 603-621. 
[19] N. Masmoudi, K. Nakanishi, From nonlinear Klein-Gordon equation to a system 

of coupled nonlinear Schrddinger equations. Math. Ann. 324 (2002), 359-389. 
[20] P. J. Pascual, S. Jimenez, L. Vazquez, Numerical Simulations of a Nonlinear 

Klein-Gordon Model, Lecture Notes in Physics 448 (1995), 211-270. 
[21] L. R. Petzold, L. O. Jay, J. Yen, Numerical solution of highly oscillatory ordinary 

differential equations, Acta Numcrica (1997), 437 483. 
[22] J. J. Sakurai, Advanced Quantum Mechanics, Addison Wesley (1967). 
[23] W. Strauss, L. Vazquez, Numerical solution of a nonlinear Klein-Gordon equa- 
tion, Journal of Computational Physics 28 (1978), 271-278. 
[24] M. Tsutsumi, Nonrelativistic approximation of nonlinear Klein- Gordon equations 

in two space dimensions. Nonlinear Anal. 8 (1984), 637-643. 



Erwan Faou, INRIA & ENS Cachan Bretagno, Avenue Robert Schumann F-35170 Bruz, 

France. • iJ-mai/ ; Erwan. Faou@inria. fr 
Katharina Schratz, INRIA & ENS Cachan Bretagne, Avenue Robert Schumann F- 

35170 Bruz, Ftance. • E-mail : Katharina. SchratzSinria.fr 



